Información sobre el estudio y objetivos

En el estudio se han seleccionado artificialmente durante 35 generaciones moscas, machos y hembras, que eran resistentes o sensibles al etanol y un grupo control con moscas seleccionadas al azar. Los grupos de moscas se seleccionaron en dos líneas diferentes. Este dataset corresponde al estudio Phenotypic and transcriptional response to selection for alcohol sensitivity in Drosophila melanogaster.

Aunque en el estudio original se analizan más factores, para esta práctica simularemos que el investigador nos ha pedido que investiguemos:

  • Main Effect Macho vs. Hembra
  • Single Effect Resistente vs Sensible
  • Interacción Sexo vs (Resistente vs Sensible)

Nota: Seguramente en este informe hemos incluido un información en exceso y alguna información quizá no fuera particularmente relevante para el investigador, sin embargo por completitud vamos a incluir un muy breve resumen de lo que haremos en cada paso y de las conclusiones que sacamos en cada uno de ellos.

Dataset

Información del dataset
Gene Expression Omnibus ID GSE7614
Organismo Drosophila Melanogaster
Número de muestras 24 (15 individuos por muestra)
Grupos Macho Control (x2 Rep. Biológicas x2 Líneas)
Hembra Control (x2 Rep. Biológicas x2 Líneas)
Macho Sensible (x2 Rep. Biológicas x2 Líneas)
Hembra Sensible (x2 Rep. Biológicas x2 Líneas)
Macho Resistente (x2 Rep. Biológicas x2 Líneas)
Hembra Resistente (x2 Rep. Biológicas x2 Líneas)
Número de características 18952
Plataforma [Drosophila_2] Affymetrix Drosophila Genome 2.0 Array
GPL1322
Tipo de muestra in situ oligonucleotide

Cargado de datos crudos

En este caso hemos descargado todos los archivos .CEL del dataset GSE7614. También hemos consultado a que grupo pertenece cada una de las 24 muestras que contiene el dataset. Hemos preparado también un archivo que contiene toda la información referente a los grupos y las muestras y lo hemos cargado en memoria.

short_name sigue la convención c(ontrol)/r(esistant)/s(ensible) m(ale)/f(emale). Los números separan las diferentes samples y las grupos de los que salieron las muestras (ver información sobre el dataset).

Luego hemos generado un ExpressionFeatureSet juntando los archivos .CEL y annotatedDataFrame que contiene el phenoData con un total 535824 filas/sondas y 24 muestras.


Control de calidad de los datos crudos

Para realizar un control de la calidad de los datos crudos vamos utilizar la libreria Array Quality Metrics y affyPLM.

Lo haremos de dos maneras diferentes, primero con los datos sin transformar y otra vez transformándolo con el logaritmo.

Sin transformar

Salida Array Quality Metrics sin transformar

Salida Array Quality Metrics sin transformar

  • En el gráfico A podemos observar un resumen de todas las muestras y el resultado obtenido del control de outliers. Vemos que según el algoritmo las muestras 1 y 17 podrían presentar algún problema por ello las tendremos en especial consideración cuando observemos los gráficos. En este caso no disponemos de la posibilidad de hacerlo, pero podríamos comprobar las imágenes de los microarrays para ver si ha ocurrido alǵun problema durante su procesado.

  • En el gráfico B observamos una agrupación por distancias de los microarrays, vemos que hay una separación clara según el sexo de la muestra, a excepción de las muestras 1 y 17 que muestran una distancia un poco atípica se encuentran más distantes del resto de muestras de machos.

  • En el gráfico C observamos los boxplots de las intensidades de cada uno de los arrays. También aparece marcado el array 17 quizá más disperso que el resto, también aparece marcado el 10, sin embargo no parece que sean muestras particularmente aberrantes.

  • En el gráfico D podemos observar un análisis por componentes principales, veremos debajo uno como el que se ha realizado durante las prácticas de la asignatura, en este caso vemos una separación clara entre los grupos macho y hembra de las muestras lo cual es esperable y deseable.

  • En el gráfico E, vemos que las muestra siguen un mismo patrón, aunque la escala hace un poco complicada esta observación, probaremos luego con una función de affysimple.

  • En el gráfico F, el valor esperable es que la mediana sea constante, sin embargo aparece una pequeña joroba en el lado derecho, pero seguramente esté relacionado con una saturación de la intensidades como indica el propio paquete.

  • En el gráfico G, podemos observar la comparación de cada uno de los arrays respecto a un pseudo-array que tiene el valor mediano a través del resto de arrays. En este caso lo esperable es que los valores se concentren a lo largo del eje M=0. Sin embargo al no estar transformados los valores aparece agrupados en el margen izquierdo haciendo difícil la visualización.

Vamos ahora a ver un PCA más completo

Vemos una agrupación clara por la condición sexo, no tanto por las condiciones experimentales.

Log-Transform

Array Quality Metrics Output

Array Quality Metrics Output

Al estar transformados con el logaritmo la forma de las gráficas cambia ya que hemos estirado la distribución. El significado de los gráficos es el mismo que en el apartado anterior.

El array 17 sigue apareciendo como un outlier, sin embargo en ninguno de los plots aparece particularmente dramático por tanto lo conservaremos.

affyPLM

Vamos ahora a utilizar la librería affyPLM vamos a ajustar un modelo a nivel de sonda y vamos a comprobar las gráficas RLE y NUSE.

En el caso del RLE vamos a encontrar una medida estandarizada de la expresión, que ha de ser similar a través de arrays, y en el caso de NUSE observaremos los residuos normalizados del modelo cuya distribución ha de ser similar a través de arrays.

No parece que ninguno de los arrays se desvíe particularmente en ninguna de las dos gráficas (quizá la c_1_m_1), si que vemos que hay una diferencia clara entre sexos. Podríamos preguntarnos en este caso si se ha cometido algún error y los machos perteneces a un batch y las hembras a otros, es decir hemos generado algún tipo de confound.

Conclusiones de calidad

  • Conservamos todos los arrays
  • Podría ser útil revisar el 17
  • Podría ser útil consultar si ha habido alguna diferencia sistemática entre los arrays de machos y hembras, batch, fecha, operario,…

Normalización

Procedemos ahora a normalizar los datos con el método RMA del paquete oligo.

Durante este paso el algoritmo:

  • Eliminará el ruido de fondo
  • Normalizará la señal entre arrays, nunca entre condiciones experimentales.
  • Colapsará todos los valores de cada grupo de sondas en un único valor por gen.

Una vez normalizado vamos a volver a realizar un análisis de calidad de los datos

Array quality metrics output después de la normalización

Array Quality Metrics Output

Array Quality Metrics Output

Los gráficos son los mismos que en apartados anteriores, ahora vemos una separación mucho más clara en el array de distancias. El array número 17 ahora es más similar al resto.

Solo existen dos detalles destacables, uno es la diferencia entre patrones de intensidad en la figura E entre las muestras machos y hembras, y la joroba que aparece en la figura F, cuando la mediana debería de ser plana, quizá sea debida a la diferencia entre sexos.

Conclusiones

Los datos parecen correctos, por tanto no es necesaria ninguna acción.

Filtrado no específico

Vamos a filtrar, lo haremos de manera manual: - Sondas de control Affymetrix - Sondas sin entrada ENTREZID - Sondas con varios genes asignados - Genes con varias sondas asignadas

Chip: Affymetrix Drosophila Genome 2.0 Array Paquete: ‘drosophila2.db’.

Quizá exista alguna manera mejor, o más óptima de tratar los casos de sondas con varios y viceversa.

De esta manera pasamos de tener 18952 features a 11294. Lo cual nos ayudará, por ejemplo, reduciendo el impacto de la correción por múltiples comparaciones.

Modelo y contrastes

Estamos ante un diseño de ANOVA 3x2, en este caso vamos a simular que estamos interesados en las comparaciones:

  • Main Effect Macho vs. Hembra
  • Single Effect Resistente vs Sensible
  • Interacción Sexo vs (Resistente vs Sensible)

Vamos a ajustar un modelo lineal con Limma y la siguiente matriz de diseño

El eje horizontal indica los términos de la matriz de diseño. El eje vertical indica el nombre de cada muestra

El eje horizontal indica los términos de la matriz de diseño. El eje vertical indica el nombre de cada muestra

Una vez ajustado el modelo calculamos los siguientes contrastes

El eje horizontal indica cada uno de los términos del modelo. El eje vertical indica cada uno de los contraste que vamos a utilizar.

El eje horizontal indica cada uno de los términos del modelo. El eje vertical indica cada uno de los contraste que vamos a utilizar.

Anotación

Anotamos los datos con el paquete drosophila2.db.

Resultados

Genes diferencialmente expresados entre condiciones

Las tablas permiten ordenar y filtrar los genes en función de diversos parámetros. Solamente se incluyen los 1000 genes más significativos por cada comparación.

Se incluyen los p-valores ajustados usando false discovery rate.

Macho vs. Hembra

Resistente vs. Sensible

Interacción Sexo y (Resistente vs. Sensible)

Volcano plot

Aquí presentamos un volcano plot con algunos extras sobre el original, con los p-valores ajustados y sin ajustar. Además se incluye una tabla con el porcentaje de genes diferencialmente expresados en cada uno de los contrastes. Se incluye el corte en p<.05 para valores ajustados y sin ajustar.

Volcano plot modificado. Las líneas horizontales indican p-valor=.05. Las líneas  verticales indican el un fold-change de 2.

Volcano plot modificado. Las líneas horizontales indican p-valor=.05. Las líneas verticales indican el un fold-change de 2.

Porcentaje de genes significativos
contrast Sin corregir p<.05 Corregido FDR p<.05
int_sex_resvssens 5.2% 0.13%
male_vs_female 86% 86%
res_vs_sens 19% 4.1%
a Redondeado a dos decimales

Diagrama de Venn

Los resultados que aparecen aquí están corregidos al p<.05, con el p-valor corregido a través de múltiples comparaciones con FDR (Bejamini Hochberg) y se han considerado los contrastes de manera separada para esta correción.

Número de genes diferencialmente expresados por contraste
res_vs_sens male_vs_female int_sex_resvssens
Down 8 5016 0
NotSig 11261 3017 11294
Up 25 3261 0
Diagrama de Venn que muestra el número de contrastes significativos encada uno de los contrastes y sus intersecciones.

Diagrama de Venn que muestra el número de contrastes significativos encada uno de los contrastes y sus intersecciones.

Tabla de anotaciones

Se ha generado (aunque no sea accesible desde el documento, si que está en el repositorio de git) con todos los genes anotados que han resultado signficativos en cada uno de los contrastes.

Archivo: ‘annotations.html’

## Warning in result_fetch(res@ptr, n = n): SQL statements must be issued with
## dbExecute() or dbSendStatement() instead of dbGetQuery() or dbSendQuery().

## Warning in result_fetch(res@ptr, n = n): SQL statements must be issued with
## dbExecute() or dbSendStatement() instead of dbGetQuery() or dbSendQuery().

Análisis de enriquecimiento

  • Hemos utilizado el paquete clusterProfiler
    • Para estos análisis hemos seleccionado todos aquellos genes cuyo p-valor corregido con FDR es menor a .1
    • Hemos realizado un análisis por contraste
    • Hemos utilizado las bases de datos KEGG y GO
    • En el caso de KEGG hemos corregido siempre con un p<.1
    • En el caso de GO hemos corregido siempre con un p<.1 y un q<.05
    • Hemos rebajado ligeramente el umbral de significación para poder mostrar los gráficos y los tipos de análisis disponibles, de otra manera esta sección quedaría vacía.
      • En un análisis real mantendríamos el umbral de signifiación fijado a priori, o marcaríamos el análisis como exploratorios.
    • Hemos realizado 10000 permutaciones en los GSE


    Enlaces a los tipos de gráficos utilizados:

En este sección utilizaremos dos tipos de análisis sobre dos bases de datos diferentes:

  • Overrepresentation y Gene Set Enrichment Analysis
    • Gene Ontology GO y Kyoto Encyclopedia of Genes and Genomes, KEGG

Interpretación de los resultados

Para los resultados de cada uno de estos análisis presentamos el mismo conjunto de tablas y resultados. Una tabla que muestra los resultados obtenidos, con información relevante sobre cada uno de ellos y un enlace a cada uno de los elementos de la GO o KEGG encontrados. Los tipos de gráficos varían según el tipo de análisis.

En el caso de KEGG cuando una ruta resulta significativa, el diagrama de la ruta se descarga y se colorea dependiendo de la expresión de los genes incluidos en la ruta para el análisis.

Gene Ontology

Over-representation analysis

Macho vs. Hembra

En este caso aparecen muchos genes relacionados con los diferentes términos de GO y la visualización se hace un poco complicada. Un upsetplot hubiese sido más útil, pero en este caso el paquete daba un error de compatibilidad que no he sido capaz de solucionar.

Resultados de sobrerrepresentación para GO

Resultados de sobrerrepresentación para GO

Resultados de sobrerrepresentación para GO

Resultados de sobrerrepresentación para GO

Resultados de sobrerrepresentación para GO

Resultados de sobrerrepresentación para GO

Resistente vs. Sensible

Resultados de sobrerrepresentación para GO

Resultados de sobrerrepresentación para GO

Resultados de sobrerrepresentación para GO

Resultados de sobrerrepresentación para GO

Resultados de sobrerrepresentación para GO

Resultados de sobrerrepresentación para GO

Interacción Sexo y (Resistente vs. Sensible)

Gene Set enrichment analysis

Macho vs. Hembra

Resultados de GSE para GO

Resultados de GSE para GO

Resultados de GSE para GO

Resultados de GSE para GO

Resistente vs. Sensible

Interacción Sexo y (Resistente vs. Sensible)

KEGG

Over-representation analysis

Contrast Rutas sobrerrepresentadas
Macho vs. Hembra 2
Resistente vs. Sensible 0
Interacción Sexo y (Resistente vs. Sensible) 4

Macho vs. Hembra

Resultados de sobrerrepresetanción para KEGG

Resultados de sobrerrepresetanción para KEGG

Resultados de sobrerrepresetanción para KEGG

Resultados de sobrerrepresetanción para KEGG

Resistente vs. Sensible

Interacción Sexo y (Resistente vs. Sensible)

Resultados de sobrerrepresetanción para KEGG

Resultados de sobrerrepresetanción para KEGG

Resultados de sobrerrepresetanción para KEGG

Resultados de sobrerrepresetanción para KEGG

Gene Set enrichment analysis

Contrast Rutas sobrerrepresentadas
Macho vs. Hembra 0
Resistente vs. Sensible 2
Interacción Sexo y (Resistente vs. Sensible) 0

Macho vs. Hembra

Resistente vs. Sensible

Resultados de GSE para KEGG

Resultados de GSE para KEGG

Resultados de GSE para KEGG

Resultados de GSE para KEGG

Interacción Sexo y (Resistente vs. Sensible)

Archivos generados

## Note: zip::zip() is deprecated, please use zip::zipr() instead
Nombre
normalized_data.csv/.Rdata
normalized_filtered_data.csv/.Rdata
top1000_male_female.csv
top1000_interaction_male_female_res_sems.csv
top1000_res_sens.csv
ora_go_interaction_male_female_res_sems.csv
ora_go_interaction_res_sens.csv
ora_go_interaction_male_female.csv
gse_go_interaction_male_female_res_sens.csv
gse_go_interaction_res_sens.csv
gse_go_interaction_male_female.csv
ora_kegg_interaction_male_female_res_sems.csv
ora_kegg_interaction_res_sens.csv
ora_kegg_interaction_male_female.csv
gse_kegg_interaction_male_female_res_sems.csv
gse_kegg_interaction_res_sens.csv
gse_kegg_interaction_male_female.csv

Referencias a los paquetes utilizados

Bolstad, Benjamin M. 2004. “Low Level Analysis of High-Density Oligonucleotide Array Data: Background, Normalization and Summarization.” PhD thesis, University of California, Berkeley.

Bolstad, Benjamin M, Francois Collin, Julia Brettschneider, Ken Simpson, Leslie Cope, Rafael A Irizarry, and Terence P Speed. 2005. “Quality Assessment of Affymetrix Genechip Data.” In Bioinformatics and Computational Biology Solutions Using R and Bioconductor, edited by R. Gentleman, V. Carey, W. Huber, R Irizarry, and S Dudoit, 33–47. New York: Springer.

Brettschneider, Julia, Francois Collin, Benjamin M Bolstad, and Terence P Speed. 2007. “Quality Assessment for Short Oligonucleotide Arrays.” Technometrics In press.

Carlson, Marc. 2016. Drosophila2.db: Affymetrix Drosophila Genome 2.0 Array Annotation Data (Chip Drosophila2).

———. 2019. Org.Dm.eg.db: Genome Wide Annotation for Fly.

Carvalho, Benilton S, and Rafael A Irizarry. 2010. “A Framework for Oligonucleotide Microarray Preprocessing.” Bioinformatics 26 (19). Oxford, UK: Oxford University Press: 2363–7. https://doi.org/10.1093/bioinformatics/btq431.

Huber, W., V. J. Carey, R. Gentleman, S. Anders, M. Carlson, B. S. Carvalho, H. C. Bravo, et al. 2015. “Orchestrating High-Throughput Genomic Analysis with Bioconductor.” Nature Methods 12 (2): 115–21. http://www.nature.com/nmeth/journal/v12/n2/full/nmeth.3252.html.

Kauffmann, Audrey, Robert Gentleman, and Wolfgang Huber. 2009. “ArrayQualityMetrics–a Bioconductor Package for Quality Assessment of Microarray Data.” Bioinformatics 25 (3): 415–6.

Luo, Weijun, Brouwer, and Cory. 2013. “Pathview: An R/Bioconductor Package for Pathway-Based Data Integration and Visualization.” Bioinformatics 29 (14): 1830–1. https://doi.org/10.1093/bioinformatics/btt285.

Ooms, Jeroen. 2020. Magick: Advanced Graphics and Image-Processing in R. https://CRAN.R-project.org/package=magick.

Pagès, Hervé, Marc Carlson, Seth Falcon, and Nianhua Li. 2019. AnnotationDbi: Manipulation of Sqlite-Based Annotations in Bioconductor.

Ritchie, Matthew E, Belinda Phipson, Di Wu, Yifang Hu, Charity W Law, Wei Shi, and Gordon K Smyth. 2015. “limma Powers Differential Expression Analyses for RNA-Sequencing and Microarray Studies.” Nucleic Acids Research 43 (7): e47. https://doi.org/10.1093/nar/gkv007.

Slowikowski, Kamil. 2020. Ggrepel: Automatically Position Non-Overlapping Text Labels with ’Ggplot2’. https://CRAN.R-project.org/package=ggrepel.

Smith, Colin A. 2019. Annaffy: Annotation Tools for Affymetrix Biological Metadata.

Warnes, Gregory R., Ben Bolker, Lodewijk Bonebakker, Robert Gentleman, Wolfgang Huber, Andy Liaw, Thomas Lumley, et al. 2020. Gplots: Various R Programming Tools for Plotting Data. https://CRAN.R-project.org/package=gplots.

Wickham, Hadley. 2016. Ggplot2: Elegant Graphics for Data Analysis. Springer-Verlag New York. https://ggplot2.tidyverse.org.

Wickham, Hadley, Mara Averick, Jennifer Bryan, Winston Chang, Lucy D’Agostino McGowan, Romain François, Garrett Grolemund, et al. 2019. “Welcome to the tidyverse.” Journal of Open Source Software 4 (43): 1686. https://doi.org/10.21105/joss.01686.

Wilke, Claus O. 2019. Cowplot: Streamlined Plot Theme and Plot Annotations for ’Ggplot2’. https://CRAN.R-project.org/package=cowplot.

Yu, Guangchuang, Li-Gen Wang, Yanyan Han, and Qing-Yu He. 2012. “ClusterProfiler: An R Package for Comparing Biological Themes Among Gene Clusters.” OMICS: A Journal of Integrative Biology 16 (5): 284–87. https://doi.org/10.1089/omi.2011.0118.

Zhu, Hao. 2019. KableExtra: Construct Complex Table with ’Kable’ and Pipe Syntax. https://CRAN.R-project.org/package=kableExtra.